{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 27,
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "x_l = 0.0\n",
    "x_r = 1.0\n",
    "n = 32\n",
    "h = (x_r - x_l) / n"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "source": [
    "def function_f(x):\n",
    "    return 1\n",
    "\n",
    "def function_u(x):\n",
    "    return -(x*x)/2\n",
    "\n",
    "def GS(A,u,Rhs):\n",
    "    v = u\n",
    "    for i in range(len(u)):\n",
    "        temp = Rhs[i]\n",
    "        for j in range(len(u)):\n",
    "            if(i != j):\n",
    "                temp -= A[i][j]*v[j]\n",
    "        v[i] = temp / A[i][i]\n",
    "    return v\n",
    "\n",
    "def L2_Error(u,v):\n",
    "    res = 0.0\n",
    "    for i in range(len(u)):\n",
    "        res += (u[i] - v[i])*(u[i] -v[i])\n",
    "    return np.sqrt(res)\n",
    "\n",
    "def interpolate(u_2h):\n",
    "    n_size = len(u_2h)\n",
    "    u_h = np.zeros((2*n_size -1))\n",
    "    for i in range(n_size-1):\n",
    "        u_h[2*i] = u_2h[i]\n",
    "        u_h[2*i+1] = (u_2h[i] + u_2h[i+1]) /2\n",
    "    return u_h\n",
    "\n",
    "def restrict(u_h):\n",
    "    u_h_size = len(u_h)\n",
    "    u_2h_size = int(u_h_size / 2) + 1\n",
    "    u_2h = np.zeros(u_2h_size)\n",
    "    u_2h[0] = u_h[0]\n",
    "    u_2h[u_2h_size-1] = u_h[u_h_size -1]\n",
    "    for i in range(u_2h_size -1):\n",
    "        u_2h[i] = (u_h[2*i-1] + 2* u_h[2 *i] + u_h[2*i+1]) / 4\n",
    "    return u_2h\n"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "source": [
    "Mat_fine_A = np.zeros((n+1,n+1))\n",
    "Rhs_fine_f = np.zeros((n+1))\n",
    "Mat_fine_A[0][0] = 1\n",
    "\n",
    "for i in range(n-1):\n",
    "    Mat_fine_A[i+1][i] = -1\n",
    "    Mat_fine_A[i+1][i+1] = 2\n",
    "    Mat_fine_A[i+1][i+2] = -1\n",
    "    Rhs_fine_f[i+1] = function_f(x_l + h*(i+1)) \n",
    "Mat_fine_A[n][n] = 1\n",
    "Mat_fine_A = Mat_fine_A /(h*h)\n",
    "Rhs_fine_f[0] = function_u(x_l)*Mat_fine_A[0][0]\n",
    "Rhs_fine_f[n] = function_u(x_r)*Mat_fine_A[n][n]\n",
    "\n",
    "###处理边界\n",
    "Rhs_fine_f[1] -= Mat_fine_A[1][0] * Rhs_fine_f[0]\n",
    "Rhs_fine_f[n-1] -= Mat_fine_A[n-1][n] * Rhs_fine_f[n]\n",
    "Mat_fine_A[1][0] = 0\n",
    "Mat_fine_A[n-1][n] = 0 "
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "Mat_coarse_A = np.zeros((int(n/2)+1,int(n/2)+1))\n",
    "Rhs_coarse_f = np.zeros((int(n/2)+1))\n",
    "Mat_coarse_A[0][0] = 1\n",
    "Rhs_coarse_f[0] = function_u(x_l)\n",
    "for i in range(int(n/2)-1):\n",
    "    Mat_coarse_A[i+1][i] = -1\n",
    "    Mat_coarse_A[i+1][i+1] = 2\n",
    "    Mat_coarse_A[i+1][i+2] = -1\n",
    "    Rhs_coarse_f[i+1] = function_f(x_l + 2*h*(i+1))*(2*h)*(2*h)\n",
    "Mat_coarse_A[int(n/2)][int(n/2)] = 1\n",
    "Rhs_coarse_f[int(n/2)] = function_u(x_r)\n",
    "\n",
    "###处理边界\n",
    "Rhs_coarse_f[1] -= Mat_coarse_A[1][0] * Rhs_coarse_f[0]\n",
    "Rhs_coarse_f[int(n/2)-1] -= Mat_coarse_A[int(n/2)-1][int(n/2)] * Rhs_coarse_f[int(n/2)]\n",
    "Mat_coarse_A[1][0] = 0\n",
    "Mat_coarse_A[int(n/2)-1][int(n/2)] = 0 "
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "print(Mat_coarse_A)\n",
    "print(Rhs_coarse_f)\n",
    "\n",
    "u_2h = np.linalg.solve(Mat_coarse_A,Rhs_coarse_f)\n",
    "u_h = interpolate(u_2h)\n",
    "\n",
    "vis = 5\n",
    "for i in range(vis):\n",
    "    u_h = GS(Mat_fine_A, u_h, Rhs_fine_f)\n",
    "\n",
    "print(L2_Error(np.dot(Mat_fine_A,u_h),Rhs_fine_f))"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "source": [
    "###进行5次迭代\n",
    "u = np.zeros((n+1))\n",
    "vis = 5\n",
    "for i in range(vis):\n",
    "    u = GS(Mat_fine_A, u, Rhs_fine_f)\n",
    "\n",
    "e_h = Rhs_fine_f-np.dot(Mat_fine_A,u)\n",
    "e_2h = restrict(e_h)\n"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "source": [
    "print(u)\n",
    "\n",
    "print(L2_Error(np.dot(Mat_fine_A,u),Rhs_fine_f))"
   ],
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "[ 0.          0.00166702  0.00279808  0.00355053  0.00404239  0.00435889\n",
      "  0.0045597   0.00468549  0.00476339  0.00481114  0.00484012  0.00485756\n",
      "  0.00486798  0.00487415  0.00487778  0.00487991  0.00488114  0.00488186\n",
      "  0.00488227  0.0048825   0.00488264  0.00488271  0.00488276  0.00488278\n",
      "  0.0048828   0.0048828   0.00488281 -0.01074219 -0.05767822 -0.13998413\n",
      " -0.25007629 -0.37454987 -0.5       ]\n",
      "0.05689678772115155\n"
     ]
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "source": [
    "Mat_coarse_A = np.zeros((int(n/2)+1,int(n/2)+1))\n",
    "Mat_coarse_A[0][0] = 1\n",
    "for i in range(int(n/2)-1):\n",
    "    Mat_coarse_A[i+1][i] = -1\n",
    "    Mat_coarse_A[i+1][i+1] = 2\n",
    "    Mat_coarse_A[i+1][i+2] = -1\n",
    "Mat_coarse_A[int(n/2)][int(n/2)] = 1\n",
    "\n",
    "r_2h = np.zeros((int(n/2)+1))\n",
    "\n",
    "r_2h = np.linalg.solve(Mat_coarse_A, 4*e_2h)\n"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "source": [
    "print(r_2h)\n",
    "print(L2_Error(e_2h, np.dot(Mat_coarse_A,r_2h)))"
   ],
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "[ 0.0004406  -0.00433815 -0.01146936 -0.02177984 -0.03568295 -0.05336533\n",
      " -0.0749051  -0.10033307 -0.12966084 -0.16289263 -0.20002991 -0.2410732\n",
      " -0.28602265 -0.33487832 -0.3250792  -0.18934929  0.        ]\n",
      "0.11209548895910082\n"
     ]
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "source": [
    "r_h = interpolate(r_2h)\n",
    "u += r_h\n",
    "for i in range(vis):\n",
    "    u = GS(Mat_fine_A, u, Rhs_fine_f)\n",
    "\n",
    "e_h = Rhs_fine_f-np.dot(Mat_fine_A,u)"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "source": [
    "print(L2_Error(np.dot(Mat_fine_A,u),Rhs_fine_f))\n"
   ],
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "0.001044989368210337\n"
     ]
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "source": [
    "for i in range(4*vis):\n",
    "    u = GS(Mat_fine_A, u, Rhs_fine_f)\n",
    "    print(L2_Error(np.dot(Mat_fine_A,u),Rhs_fine_f))\n",
    "print(np.dot(Mat_fine_A,u)-Rhs_fine_f)"
   ],
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "0.00024336891930645596\n",
      "0.00022379251701409283\n",
      "0.00020653198398767937\n",
      "0.00019122683742338927\n",
      "0.00017758657352207768\n",
      "0.00016537443117220211\n",
      "0.00015439551829612316\n",
      "0.00014448797067896974\n",
      "0.00013551626652162784\n",
      "0.00012736610319555067\n",
      "0.00011994042575788242\n",
      "0.00011315631548509664\n",
      "0.00010694252792579144\n",
      "0.00010123752643521264\n",
      "9.598789717184268e-05\n",
      "9.11470605852044e-05\n",
      "8.667421590519695e-05\n",
      "8.253347110261093e-05\n",
      "7.869312250236329e-05\n",
      "7.512505664446867e-05\n",
      "[ 0.00000000e+00  2.53904681e-07  1.70896638e-07 -1.31373422e-07\n",
      " -7.06855528e-07 -1.59218614e-06 -2.79289922e-06 -4.26991015e-06\n",
      " -5.92970778e-06 -7.62169002e-06 -9.14517816e-06 -1.02669627e-05\n",
      " -1.07481102e-05 -1.03766777e-05 -9.00142472e-06 -6.56095968e-06\n",
      " -3.10316280e-06  1.20886501e-06  6.10596632e-06  1.12364342e-05\n",
      "  1.62008462e-05  2.05935218e-05  2.40450332e-05  2.62602475e-05\n",
      "  2.70472633e-05  2.63341071e-05  2.41718936e-05  2.07249928e-05\n",
      "  1.62503324e-05  1.10690842e-05  5.53454211e-06  0.00000000e+00\n",
      "  0.00000000e+00]\n"
     ]
    }
   ],
   "metadata": {}
  }
 ],
 "metadata": {
  "orig_nbformat": 4,
  "language_info": {
   "name": "python",
   "version": "3.7.10",
   "mimetype": "text/x-python",
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "pygments_lexer": "ipython3",
   "nbconvert_exporter": "python",
   "file_extension": ".py"
  },
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3.7.10 64-bit ('base': conda)"
  },
  "interpreter": {
   "hash": "2c21737055494eda6dd7820cee4b5069ee6eef9e680dd0a23b71fc1db8cdc626"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}